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Abstract 

The paper develops high-order accurate physical-constraints-preserving finite dif¬ 
ference WENO schemes for special relativistic hydrodynamical (RHD) equations, 
built on the local Lax-Friedrich splitting, the WENO reconstruction, the physical- 
constraints-preserving flux limiter, and the high-order strong stability preserving 
time discretization. They are extensions of the positivity-preserving finite difference 
WENO schemes for the non-relativistic Euler equations [21], However, developing 
physical-constraints-preserving methods for the RHD system becomes much more 
difficult than the non-relativistic case because of the strongly coupling between the 
RHD equations, no explicit formulas of the primitive variables and the flux vectors 
with respect to the conservative vector, and one more physical constraint for the 
fluid velocity in addition to the positivity of the rest-mass density and the pressure. 
The key is to prove the convexity and other properties of the admissible state set 
and discover a concave function with respect to the conservative vector instead of 
the pressure which is an important ingredient to enforce the positivity-preserving 
property for the non-relativistic case. 

Several one- and two-dimensional numerical examples are used to demonstrate ac¬ 
curacy, robustness, and effectiveness of the proposed physical-constraints-preserving 
schemes in solving RHD problems with large Lorentz factor, or strong discontinu¬ 
ities, or low rest-mass density or pressure etc. 

Key words: Finite difference scheme; Physical-constraints-preserving; High-order 
accuracy; Weighted essentially non-oscillatory; Relativistic hydrodynamics; 

Lorentz factor. 
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1 Introduction 


The paper is concerned with developing high-order accurate numerical methods for spe¬ 
cial relativistic hydrodynamical (RHD) equations. In the laboratory frame, the {d -|- 1)- 
dimensional space-time RHD equations may be written into a system of conservation laws 
as follows 


dU 


i=l 


dxi 


( 1 . 1 ) 


where U and are the conservative vector and the flux in the Xj-direction, respectively, 
dehned by 


U =(^D,mi, ■■■ , 

Fi =[Dvi,miVi+■ ■ ,mdVi+p5d,i,rn^ , i = l,---,d, 

with the mass density D = pW, the momentum density vector m = DhWv, and the 
energy density E = DhW — p, respectively. Here p, p, and n = (ui, • • • , Vd)'^ denote the 
rest-mass density, the kinetic pressure, and the fluid velocity respectively, W = l/\/l — 
is the Lorentz factor with v = (ni + - • ■ + VdY^'^, and h denotes the specihc enthalpy dehned 
by 

h = 1 -|- e H—, (1-2) 

P 

with units in which the speed of light is equal to one, and e is the specihc internal energy. 

The system fll.ip has taken into account the relativistic description of huid dynamics 
where the huid how is at nearly speed of light in vacuum and appears in investigating 
numerous astrophysical phenomena from stellar to galactic scales, e.g. formation of black 
holes, coalescing neutron stars, core collapse super-novae. X-ray binaries, active galactic 
nuclei, super-luminal jets and gamma-ray bursts, etc. However, it still involves highly 
nonlinear equations due to the Lorentz factor so that its analytic treatment is extremely 
difficult. A powerful and primary approach to improve our understanding of the physical 
mechanisms in RHDs is through numerical simulations. Comparing to the non-relativistic 
case, the numerical difficulties are coming from strongly nonlinear coupling between the 
RHD equations, which leads to no explicit expression of the primitive variables V = 
(p, and the hux vector F, in terms of U, and some physical constraints such as 

p > 0, p > 0, E > D, and 1 > n etc. Its numerical study did not attract considerable 
attention until 1990s. 

* Corresponding author. 
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The first attempt to numerically solve the RHD equations was made by using a hnite 
difference method with the artificial viscosity technique in Lagrangian or Eulerian co¬ 
ordinates, see [?T]l32ll48lim] . After that, various modern shock-capturing methods were 
gradually developed for the special RHDs since 1990s, e.g. the HLL (Harten-Lax-van 
Leer-Einfeldt) method |13], the two-shock approximation solvers [T1I8] . the flux corrected 
transport method HU, the Roe solver [13], the HLLC (Harten-Lax-van Leer-Contact) 
approximate Riemann solver | 37 |, the flux-splitting method based on the spectral decom¬ 
position [10], Steger-Warming flux vector splitting method [69], and the kinetic schemes 
[Fni271139] ■ Besides those, high-order accurate schemes for the RHD system were also 
studied, e.g. ENO (essentially non-oscillatory) and weighted ENO methods [9)1601146] . the 
piecewise parabolic methods [331)38] , the space-time conservation element and solution el¬ 
ement method [iQ], and the discontinuous Galerkin (DG) method [12], the Runge-Kutta 
DG methods with WENO (weighted ENO) limiter [70], the direct Eulerian GRP schemes 


j , the adaptive moving mesh methods 
hnite volume local evolution Galerkin method 
early review articles 



, and genuinely multi-dimensional 
The readers are also referred to the 


The above existing works do not preserve the positivity of the rest-mass density and the 
pressure and the bounds of the huid velocity. Although they have been used to simulate 
some RHD hows successfully, there exists the big risk of failure when they are applied to 
RHD problems with large Lorentz factor, or low density or pressure, or strong disconti¬ 
nuity, because the negative density or pressure, or the larger velocity than the speed of 
light may be obtained so that the calculated eigenvalues of the Jacobian matrix become 
imaginary, in other words, the discrete problem becomes ill-posed. In practice, the non¬ 
physical numerical solutions are usually simply replaced with a “close” and “physical” 
one by performing recalculation with more dihusive schemes and smaller CEL number 
until the numerical solutions become physical, see e.g. [6D[22] . Obviously, such approch is 
not scientihcally reasonable to a certain extent, and it is of great signihcance to develop 
high-order accurate numerical schemes, whose solutions satisfy the intrinsic physical con¬ 
straints. Recently, there exist some works on the maximum-principle-satisfying schemes 
for scalar hyperbolic conservation law [62167], the positivity-preserving schemes for the 
non-relativistic Euler equations with or without source terms [66)1631164] . the positivity¬ 
preserving well-balanced schemes for the shallow water equations [S3], the positivity pre¬ 
serving semi-Lagrangian DG method for the Vlasov-Poisson system [TT], and Lagrangian 
method with positivity-preserving limiter for multi-material compressible flow [5]. A class 
of the parametrized maximum principle preserving and positivity-preserving flux lim¬ 
iters were also well-developed via decoupling some linear or nonlinear constraints for 
the high order accurate schemes for scalar hyperbolic conservation laws [59)129] , nonlinear 
convection-dominated diffusion equations [25], compressible Euler equations [53] and ideal 
magnetohydrodynamical equations [6] as well as simulating incompressible flows [55] • A 
survey of the maximum-principle-satisfying or positivity-preserving high-order schemes is 
presented in [65] . 


The aim of the paper is to do the first attempt in the aspect of developing the high-order 
accurate physical-constraints-preserving finite difference schemes for special RHD equa- 
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tions (inj. Such attempt is nontrivial in comparison of the non-relativistic case, because 
of the strongly nonlinear coupling between the RHD equations due to the Lorentz factor, 
no explicit formulas of V and Fi with respect to [/, and one more physical constraint for 
the fluid velocity in addition to the positivity of the rest-mass density and the pressure. 
The key will be to prove the convexity and other properties of the admissible state set 
and discover a concave function with respect to the conservative vector U instead of the 
pressure, which is an important ingredient to enforce the positivity-preserving property 
for the non-relativistic case. The paper is organized as follows. Section [2] discusses the 
admissible state set and its properties of the special RHD equations. They play a piv¬ 
otal role in studying the physical-constraints-preserving property of numerical schemes. 
Section |3] presents the high-order accurate physical-constraints-preserving hnite difference 
WENO schemes for the RHD equations. Section ITT] considers detailedly one-dimensional 
case with spatial discretization in Section 13.1.11 and time discretization in Section 13.1.21 
Section 13.21 gives the 2D extension of the above scheme and apply it to the 2D axisym- 
metric case. Section 0] gives several ID and 2D numerical experiments to demonstrate 
accuracy, robustness and effectiveness of the proposed schemes for relativistic problems 
with large Lorentz factor, or strong discontinuities, or low rest-mass density or pressure 
etc. Section |5] concludes the paper with several remarks. 


2 Properties of RHD equations 


This section discusses the admissible state set and its properties for the RHD equations 
fll.ip . Throughout the paper, the equation of state (EOS) will be restricted to the T-law 

p={r-l)pe, (2.1) 

with the adiabatic index T G (1,2]. Such restriction on T is reasonable under the com¬ 
pressibility assumptions, see [7]. 

The RHD system fll.ip with fl2.ip is identical to the d-dimensional non-relativistic Euler 
equations in the formal structure, and also satishes the rotational invariance and the 
homogeneity as well as the hyperbolicity in time, see [69]. The momentum equations in 
fll.ip are only with a Lorentz-contracted momentum density replacing pvi in the non- 
relativistic Euler equations. When the fluid velocity is much smaller than the speed of 
light (i.e. n 1) and the velocity of the internal (microscopic) motion of the fluid particles 
is small, the system fll.ip reduces to the non-relativistic Euler equations. The {d + 2) real 
eigenvalues of the Jacobian matrix Ai(U) = dFi(U)/dU for fll.ll) with fl2.ip are 

(1) _ Vi{l - cl) - CsW-^^Jl-vf - {v^-vf)cl 

* “ 1-v^cl 

\( 2 ) _ _ _ 

Aj^ — — '^i — D; 

, (d+ 2 ) _ d(1 - cl) + CsW-^^l-vf - {v^-vf)cl 
' “ 1-v^cl 
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for i = I,-- - ,(i, where the local sound speed However, relations between 

the laboratory quantities {D, rrii, and E) and the quantities in the local rest frame (p, 
Vi, and e) introduce a strong coupling between the hydrodynamic equations and pose 
more additional numerical difficulties than the non-relativistic case. For example, the flux 
vectors Fi and the primitive variable vector V := {p,v,p)'^ can not be formulated in 
explicit forms of the conservative vector U, and some constraints (e.g. p > 0, p > 0, and 
1 > n, etc.) should be fulhlled by the physical solution U. 


Definition 2.1 The set of (physically) admissible states of the RHD equations (II.Ih with 
(12.ip is defined by 


g = = {D,m, £)’■[ p(U) > 0, p(U) > 0, 1 > ti(t/)} . 


( 2 . 2 ) 


Such dehnition is very natural and intuitive but unpractical when giving the value of the 
conservative vector U, because there is no explicit expression of the primitive variable V = 
(p, n,p)^ in terms to U for the system (11.ip . It is this feature that makes the discussions 
on the admissible state and the physical-constraints-preserving schemes presented later 
nontrivial or even challenging for RHD equations (11.11) . In practical computations, for 
given U = {D, m, Efi", one has to (iteratively) solve a nonlinear algebraic equation such 
as 

E + p = DW+,^^pW^, (2.3) 

by any standard root-hnding algorithm to get the pressure p{U). Then Vi and p are 
sequentially calculated by 


= sfipiur " Dp-vHU). (2,4) 

Note that the Lorentz factor W in (12.3p has been rewritten into (1 — m?/[E + pfi) 
with m := {m\ + • • • + . Existence of the unique positive solution to the pressure 

equation (12.3p is given in the proof of Lemma 12.11 It is worth mentioning that different 
from the non-relativistic case, such positive solution p(J7) is not a concave function of U 
generally, see Fig. 12.11 


A practical and equivalent dehnition of Q is given as follows. 

Lemma 2.1 The admissible set Q defined in (12.2p is equivalent to the following set 


= |t/ = {D, m, EY' D > 0, q{U) := E — V + m? > 0 1. 


(2.5) 


Proof (i) : Prove that U E Qi when U E Q. When U = {D,m, EY satisfy the 
constraints p{U) > 0, p{U) > 0, and v(U) < 1, it is not difficult to show 


D = 


P 




>0, E = 


ph 


1 -n2([/) 


p> ph-p^ p(l -P e) > 0, 
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Fig. 2.1. The function <^(A) := + (1 — A)?/®) — Xp{U^^^) — (1 — A)p(Ll®) with 

= (2,1.2, S)"*" e G and = (2,5,35)'^ G Q. The value of ip{X) is always less than 
zero when A G (0,1). 

and 


{d^ + 


m = 


ph 


1 — 


P 




ph 


— n 
1 

1 — 

V<1 1 

> 


2 . +p^ - 2p—^ 


1 — 


P 


phv 
1 — 

/ phv 


{ph — p) — p^ — p^v 


1 — 1 — \1 — v"^ 

OJ 1 


1 — 


1 — 

p2 (1 + e)2 _ ^2 _ p2] Ip 1^2 + er(2 - r)) ' 0 


2 /I I 2 2 ^ 

p (1 + e) - p - p V 
re(i, 2 ] 


Thus q{U) = E — \/D‘^ + m? > 0 so that U G G\. 

(ii) : Prove that U E G when U E Gi- Consider the function of p dehned by 


$(p) ;= 


m 


+ D, 


E + p \ {E + pY T-I 


™ +J^-E, pelO.+oo), 


with U satisfying that D > Q and E > Obviously, $(p) G C^[0, +oo), and 


$'(p) = 


m 


D 


> 1 


m 


>0, \/p E [0, +oo), 


r-1 {E + pY\^ ^{E + pY-m? j ^ {E + pY 

when E > YD'^ + m? and T G (1,2]. Thus $(p) is a strictly monotonically increasing 
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function of p in the interval [0, +cxd). On the other hand, one has 



and lini <h(p) = +oo because 

p^+oo ' ' 



> 0 . 


Thanks to the intermediate value theorem and the monotonicity of $(p), there exists a 
unique positive solution to the equation $(p) = 0, which is equivalent to the equation 
fl2.3p . Denote this positive solution hy p{U). Substituting this into fl2.4p may give 



by using the conditions that D > 0 and E > \/Thus U E Q and the proof is 


I 


completed. 


Remark 2.1 Comparing to Q defined in fl2.2p . the constraints on conservative variables 
in the set Qi is much easier to he verified when the value of U is given. 

With the help of equivalence of the admissible state set Q in Lemma 12.11 we can further 
prove that it is a convex set. 

Lemma 2.2 The admissible set Qi is a convex set. 

Proof To show that the set Qi is convex, one has to prove that for all A in the interval 
[0,1], and all and in the set the 

point + (1 “ A)t/® =: G Qi also belongs to Q\. 

Because G Gi, one has 


L)(^) = AD(^) + (1 - A)DW > 0, 


and 


+ (1 - 




2 
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Here the Minkowski inequality for both vectors Xm^i \ • • • , and ((1—, (1— 

X)mf\--- ,(1 —and the triangle inequality lAm-^^ + (1 — A)m'°^ | < A|m'^^| + (1 — 
A)|m'°^| have been used respectively. Thus AC/d) + (1 — A)L/'® G for all A G [0,1]. The 
proof is completed. | 


Remark 2.2 The proof of Lemma \2.2\ implies that the function q{U) defined in fl2.5p is 
concave. Moreover, the function q{U) is also Lipschitz continuous with respect to U and 
satisfies 


Uc/d))-g(C/W)| < Ed)_^(0) + + (md))^ - + (m(0))^ 


< 


< 


£;( i ) _ ^{ o )| + (^( 1 ) _ ^{ o ))2 ^ ^ _ ^( 0)^2 

\ i=l 


r 

= V2 


(T)(i) - T)(o))^ + ^ 

f/d) - (7(0) 


2 = 1 


( 2 . 6 ) 


for any f/(o) = (T)(o), m(o), ^(o))^ G M“'+2 jy(i) = (£>(1)^ ^(i)^ ^(i))T ^ ^d+ 2 ^ 

the inequality a + h < \j2{o? + b‘^) has been used. The concavity and Lipschitz continuity 
of q{U) will play a pivotal role in designing our physical-constraints-preserving schemes 
for the RHD equations fll.ip . 


By means of the convexity of Q, the following properties of Q can further be verihed. 


Lemma 2.3 Assume U E Qi, then 


(i) Xu G Qi, for all A > 0. 

(a) TU G Qi, where T = diag{l,Tdyid^ 1} and T^xd is the dx d rotational matrix. 

(Hi) U ± a~^Fi(U) G Qi for all real number a > Qi, i = 1, ■ ■ ■ ,d, where Qi is the spectral 
radius of the Jacobian matrix Ai(U), i.e. 

luil (1 - cl) + CsW~^Jl - vj - - v})cl 

Qi ■ 2~2 ■ 

Proof The verihcation of hrst two properties are omitted here because they can be di¬ 
rectly and easily verihed via the dehnition of Qi. 

The following will prove the third conclusion [Hi) that if [/ G Qi, then = 

(7^ ■.= U E a ^FifU) G Qi. It is nontrivial and requires several techniques thanks to no 
explicit formulas of the hux Fi{U) in terms of U. 


For all U E Qi and the perfact gas fl2.ll) with T G (1, 2], one has 


Tp 




Tp 


ph P+ friP 


< T - 1 < 1, 


(2.7) 
































so that the inequalities 


and 


W ^— vf — — vf)c^ > W ^\/l — = 1 — 


> 


hold. They imply 


1 - n2c2 




kiKi-t-yt 


W ‘^Cs (1 - |ni| Cs) 


^1 


> 


|ni|(l - cl) + Cs(l - |ni| (1 - cl) + c^(l - 

hh-^C, (1 - |ni|c,) TI/-2, 


IT-"c, ^ lT-"c, 


|ni| (1 - c2) + 0^(1 - nf) luil+c^ l + c, 


so that one has 


Thus one gets 


1±^>1 

a 


M>i_hl>Td£i>o, 


a 


Qi 1 + Cs 


( 2 . 8 ) 


Vl 


= D (^1 ± ^ j > 0, 

7-± 7- I "^1 t LTI72 N , LTI72^1 ® PhCs E3 

E = E ±.— = iphW — p) E phW — > -- p = p 

a a 1 + Co 


Cs(l + cj 


iMli 

> p 


T-l + x^^^ 


- 1 =p 


1 - VT^ \ rG(i,2] 

' > 0 , 


T - 1 + x/T^ 


and 


{D^y + f - {E^y 

= fl ± (d^ + ± ^ (m, - Ev,) (l ± ^2.) + ?!(! _ 

= ('l±—) W'^ \p^ + p'^v'^ - {ph - p)"^] + p'^ (l ±—'] -p^ + ^ 
V a / L J V a / 



5 ) 

= (l± 

-) 

V 

a J 

= (l± 

-) 

V 

a J 

A 

1 

1^11 


^>1 

here f|2.8p and p" 



p +p - p + 


p 


T - 1 


+ p 


- 1 




P +P - {P + 


p 


T - 1 




(2,9) 


to the following quadratic equation 

(1 - v^cDgj - 2 |ni| (1 - + n?(l - c^) - c^(l - v^) = 0, 

which is equivalent to 

(l - P?) cl = IT^ (pi - Iml)^ (1 - cl). 


( 2 . 10 ) 
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It implies that qi < 1. With the help of fl2.9p and (12.1 Op . one has 



±^2 


- 


P 


p 


(1 - (r -1) 

1 


_ 

r-1 (r-i)"_ 

r -1 - c? 


+ p [- 

\Qi 


p 


(1 - (r -1) 

p2 1 - 2 r rG(l,2] 

< 0 . 


' 1 2p' 

r-1 p 


(i-c=)(r-i) h 


Thus U ± a~^Fi{U) G Oi- Combining the above deduction and the property (ii) may 
verify U ± a~^Fi(U) E ior i = 2, ■■■, d. For example, in the case of d = 3, T may be 
taken as 


•— 


1 0 
0 cos 0 sin cj) 
0 — sin d 


0 0 0 
sin 6 sin (j) cos 0 0 
cos 6 0 0, 


0 — cos 9 cos 0 — sin 6 cos 0 sin 0 0 


0 0 0 0 1 


for all 9 E [0, 27 r), 0 E [0, vr]. Using the property (ii) that Tq^^U E Qi and the above proof 
of the property (Hi) gives 


Te,^U±a-^Fi{Te,^U) E Qi, 

where a is not less than qi. Since is also a rotational matrix, it holds that 

± a-^FQTg^^U)) =U± a-^T^^^^FQTg^^U) E Q,. 

With the help of the rotational invariance of the system (11.11) . see [69], one has U ± 
a~^Fi(U) E Qi for Va > Qi, i = 2, 3. The proof is completed. | 

It is worth emphasizing that Lemma 12.31 also plays a pivotal role in seeking the physical- 
constraints-preserving schemes. 


3 Numerical schemes 


This section gives the physical-constraints-preserving hnite difference WENO schemes for 
the RHD system (II. ip with the EOS (12.ip . 
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3.1 One-dimensional case 


This subsection first discusses numerical discretization of the ID RHD equations in the 
laboratory frame 


dU dFi{U) 
dt ^ dx 


(3.1) 


where 

U = (D, mi, E)'^, Fi = {Dvi, miVi + p, mi)^. 


3.1.1 Spatial discretization 

Let us divide the space into cells of size Ax, and denote the jth cell by Ij = (xj_i,Xjj^i^, 
where = ^{xj + Xj+i) and Xj = jAx, j E Z. 

A semi-discrete, (2r — l)th-order accurate, conservative hnite difference scheme of the ID 
RHD equations fl3.ip may be written as 

(3-2) 

where Uj{t) ~ U{xj,t) and the numerical flux Fj_^_i is consistent with the flux vector 
Fi{U) and satisfies 

. - -Fj-.) = 9.Fi(C/)U, + 

Definition 3.1 The scheme fl3.2p isphysical-constraints-preservingifUj{t)+AtC{U{t)-,j) G 
Q for all j G Z, under a CFL-type condition for At when Uj{t) G Q for all j. 

The third property of Lemma 12.31 has implied that there at least exists a physical- 
constraints-preserving scheme for the ID RHD system fl3.ll) . An example is the hrst-order 
(i.e. r = 1) accurate local Lax-Friedrichs scheme with the numerical flux 

Fj+l = -PflJ := i(Fi(t/j) + F.(t/,+.) - (t/,+, - U,)), 

where the viscosity coefficient cKj+i := max {pi ([/j_r+i), • • • ,gi{Uj+r)}- In practical 
computations, the above Q(j+i niay te replaced with 

Oj+l = 'll max {pf+i®, Pi {U j-r+l) , • • • , Pi ([/ j+r)} • 

see [2], where the parameter i? is typically in the range of 1.1 to 1.3 and controls the 
amount of dissipation in the numerical schemes while is the spectral radius of the 

Roe matrix Ai(t/j, l/j+i) approximating the Jacobian matrix Ai[U), see [13] . 
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The physical-constraints-preserving property of the above local Lax-Friedrichs scheme is 
shown below. 


Lemma 3.1 IfUj G G for all j, then under the CFL-type condition 


At < 


Ax 


2 maxa,-, 1 ’ 

j 3+2 


(3.3) 


one has 


U 


±,LLF 


:= U, T 


2 At - LLF ^ 


and Uj{t) + Ate {U{t)]j) = | G G for all j. 

Proof Note that can be rewritten as 

= U, T + Pi(t/,±i) T (t/jii - U 


= 1 


At\ 


At 
2 Ax 


UiT 





Uj±i T 


) 


which is a convex combination under the CFL-type condition fl3.3p . Utilizing the property 
(iii) of Lemma [2.31 implying that 


J, ^ -flit/,) Fi(V±i) , 

Uj T ) etj±l ^ ^ ^ Gy 






and the convexity of the set G in Lemma 12.21 may complete the proof. 


I 


The remaining task is to develop higher-order (i.e. r > 1) accurate physical-constraints- 
preserving finite difference scheme for the ID RHD equations fl3.ip . To finish such task, the 
idea of the positivity-preserving hnite difference WENO schemes for compressible Euler 
equations in [66] and the flux-limiter in [21] is borrowed here. For the sake of convenience, 
the independent variable t will be temporarily omitted. 


Given point values {Dj}, for each j, calculate 

lit ■■= ^ [Uk ± , J - r + 1 < < J + r, 

which may be considered as the point values of both local Lax-Friedrichs type splitting 
functions 

]^{u{x) ±a-l,_FGU{x))y 

If dehne the functions by 

1 / \ ][ 
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then become the cell average values of i(x) over the cell Ik because 


Ax 


r^k+^ , 

/ Hf 

^+2 


j — r + \<k<j + r. 


Based on these cell-average values, using the WENO reconstruction [2^IHi] may get high- 

order accurate left- and right-limited approximations of .{x) at the cell boundary 

2 


X 




, denoted by H 


,WENO 


, „-,WENO 

and 


respectively. After then, the numerical flux of a 


(2r — l)th-order accurate hnite difference WENO scheme of the ID RHD equations fl3.ip 
is 

/ , TTTT-yAT^ ,TT-r-,^T^\ 

(3.4) 


- -WENO 


rj-+,WENO rr-,WENO \ 

^j+lL - 


J 


Practically, for the system of conservation laws, a better and stable way to derive those 
left- and right-limited approximations is to impose the WENO reconstruction on the 
characteristic variables by means of their cell-average values 


wt := R~llH 


k 5 




j — r + l<k<j + r — 1, 
j — r + 2<k<j + r, 


where Rj+i is the right eigenvector matrix of the Roe matrix Ai{Uj, Uj^i). After having 
the left- (resp. right-) limited WENO values of the characteristic variables at the cell 
boundary x^_^_i, denoted by (resp. then calculate 


„-+,WENO „ ti;-+:WENO 


rr-,WENO B tj^-,WENO 

H. I „ = R,:iW 1 „ . 

3+n,R J+2 3+n,R 


Generally, the high-order accurate hnite difference WENO schemes fl3.2p with the numer¬ 
ical hux given in (13.4p is not physical-constraints-preserving, that is to 

say, it is possible to meet Uj{t) + AtC {U{t)]j) = \ ^jy+.WENO jy-,WENO^ ^ where 

jy±,WENO ._ ^ Thus for some demanding extreme problems, such as their 

solutions involving low density or pressure, or very large velocity or the ultra-relativistic 
how, these high-order schemes always easily break down after some time steps due to 
the nonphysical numerical solutions {Uj ^ Q). To cure such difficulties, the positivity¬ 
preserving hux limiter [2T] for non-relativistic Euler equations may be borrowed and ex¬ 
tended to our RHD case. Because of the dehnition of Qi in Lemma 12.11 and the properties 
of q{U) shown in Remark 12.21 the resulting physical-constraints-preserving hux limiter 
may be formed into two steps as follows in order to preserve the positivity of D{U) and 
q{U). The hux limiter such as the parametrized hux limiter [5^29ll54j can also be extended 
to our RHD case in a similar way. 

Before that, two sufficiently small positive numbers Ed and Eq are hrst introduced (taken 
as 10“^^ in numerical computations) such that > 0 and g([/^’^^^) > > 0 

for all j. It is true because Lemma [SH] tells us that the mass-density > 0 and 

^(j^±,LLF) ^ g where denotes the hrst component of 
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Step I: 

-WENO 

‘"iH 


Enforce the positivity of D{U). For each j, correct the numerical flux 
as 






(3.5) 


where denotes the lih. component of Fj+:^ and 9o,j+i/2 = j+ 1/21 j+ 1 / 2 } 

with 




( n±:LLF N u n±.LLF ni.WENOx -r n±.WENO ^ 

), <^z), 


j+iT 




J + 5 =F -2 

otherwise. 


Step II: Enforce the positivity of g([/). For each j, limit the numerical flux 

as 




:=(i 


^Q.i+5 


, - LLF „ - D 


where 1 

HU~r 2 


mm <,9^.^1,9 .^i 
q,j+2 9J+2 


and 


(3.6) 


e 


± 



if 

otherwise. 


It is worth emphasizing that the above limitting procedure is slightly different from that 
in [2], because only the first component of the numerical flux vector is detected and 
limited in Step I. Moreover, the finite difference schemes fl3.2p with the numerical flux 

..-V ..-V PCP 

= Fj^i given in (13.bh is obviously consistent with the ID RHD equations (13.ip 
and physical-constraints-preserving (see Theorem 13.ip . and maintains (2r — l)th-order 
accuracy in the smooth region without vacuum (see Theorem I3.2p . 

Theorem 3.1 Under the assumption of Lemma \3.1{ if D ±,LLF ^ Q J^±,LLF^ ^ Q 

all j, then e G, and Uj{t) + AtC {U{t)-j) = | G G 

for allj, where := Uj =F 

Proof According to Lemma 13.11 and the previous flux limiting procedure, one has 0 < 
j+i, < 1, and there exist two sufficiently small positive numbers Ed and Eq such 

that > 0 and > Eg > 0 for all j. 

Substituting fl3.5p into fl3.6p gives 
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where 0 < 9j^i = According to the dehnition of 

inequality 


(1 - 


,WENO 


> ^D, 


always holds. Thus one has 






-,LLF , n n+.WENO 

-,+ i +^. + i^,+l 


^^■+1 f n Q+ ^ n+TLF , n+,WENO^ , / i 


> 


^Ai+2 




D 


+,LLF 


9t 


'£d + I 1 


D,j+i 




7~\-1-,LLF V. . ^ 

DJ' > Ed > 0. 


On the other hand, similarly, making use of the concavity of q(U) gives 
9 (c/y'i"'’) = 9 ((1 - 

> (1 - (u;^Y) + ««+i9 (c/J;?) 


> 


n+ ^9 


qJ +2 
9. 




<id+i 


,D 

q+h 


9 


+ 1 - 


qd+2 

C+r 


9 (t'kr 


Sa+ 1 — 


q^i+h 


9 I 


>£g>0. 


With the equivalent definition of the admissible state set Q in Lemma [2Tl one knows that 


E Q. Similarly, one can also prove that Uj G Q. The proof is completed. | 


-,PCP 


The following is to check the accuracy of the schemes (13.2p with the numerical flux .F _,_i = 


-PGP 


.^,+1 given in fl3.6p . In fact, the above flux limiting procedure implies that 


-PGP 


WENO 




< (1 - «,+* 


-LLF -WENO 


Thus if 


l-9^.^^=0(Ax^’^-^), (3.7) 

then the schemes (13.21) with the numerical flux = ^j+i is (2r — l)th-order accurate 
because both and qj-q bounded in smooth regions. 

Theorem 3.2 Assuming that the exact solution U{x) is smooth and satisfies that D{x) > 
£ > 0 and q(U{x)) > e > 0 for all x, where e > -A:;^max{eD,Sq}, and the approximate 
solution Uj E Q and Uj = U{xj) + 0{Ax‘^'^~^) for all j and sufficiently small Ax, then 

0%,^ = 1 + 0{Ax^^-^), = 1 + 0{Ax^^-^), 
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and fl3.7p hold for the given (2r — l)th-order accurate numerical flux under the 

CFL-type condition 

, wAx 

At < --, 3.8 

/ maxa,-, 1 

j 3+2 

where w is any positive constant less than one. 

Proof Only estimations of i and .i are given below, while other cases are similar 
and will be omitted here. 


Before that, some basic conclusions are hrst listed as follows. 
• Lemma [3.11 implies that 


2At -LLF ^ 


(3,9) 


under the CFL-type condition 
Since 




-WENO 




= 0{Ax), 

where the vector function H{x) is implicitly dehned by 

1 


= 0{Ax 


2r-l^ 


one has 


C/^ ; = C/+'LLF ^ 


2At 


Ax 

^ ,y-+,LLF 2At 

^ Ax 

2 At - WENO 

= U j T —F ,.1 

3 Ax ^ + 2 




WENO 


^3H 


+ OiAx 


2r-l^ 


C>(Ax"’-^) = [/+’WENO ^ 0(^^^2r-U 


It further yields 


jje _ ^-l-,LLF 


J 3 


< 


- u, 


+,LLF 


g(C/0 - g(C/+’^^^) <V2 


-,WENO 


q{U^) - g([/; 


< 


,WENO^ 


U--U] 

< V2 


u,-u] 

,WENO 


= C>(Ax), 

LLF 


= C>(Ax). 


= 0{Ax 


2r--l3 




,WENO 


= C>(Ax""-'), 


(3.10) 

(3.11) 

(3.12) 

(3.13) 


where (12. 6 p has been used. 


(i) : Estimate 0); ., i. The case of that 9t, ., i = 1 is trival. Assuming 0 < 6*); ., i <1, 
D,] + ^ D,3 + j ^ — D,J + ^ ’ 
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which implies < Sd, then 


- Ed 


Ed — Dj 


1 _ a+ _ 1_-_ _ _ 

D,i+i ^+,LLF _ ^+,WENO ^+,LLF _ jj-i 


,WENO 


D 


j 


DJ 


,WENO 


< 


Ed — D 


H,WENO 


3 


j 


- Ed 


Since Dj is a (2r — l)th-order accnrate approximation to D{xj), one has Dj = D{xj) + 
0{Ax‘^'^~P > e + 0{Ax‘^''~p > Thanks to fl3.9l) . one gets 

-Ed = {1- w)D^ + w {d^ - {^i+i} J -eD>{l- wf- - > 0, 

from zero. Then, one only needs to show 


I LLP 

which shows that Dj’ — is bonnded away f 

r-)+,WENO ^ 0{Ax‘^'^~p. Note that fl3.10p implies 


Ed — D 


Dj - Ed = Dp^^^ - eD + 0{Ax) >{1 -w)^-ed + 0{Ax) - w)^- ed^ > 0. 


Thns 

n+.WENO 

Ed - Dj 

where fl3.12p has been nsed 


= Ed- D+’Weno < 


,WENO ^ 


(ii) : Estimate 6^., i. Similarly, only consider the nontrival case that 0 < i < 1 

9J+2 J — 

, . ._L WT?MOk _ 


HiJ~ 2 

implying q{Up^^^^) < Eg. Thns 


Because Uj is a (2r — l)th-order accurate approximation to U{xj), q{Uj) is also a (2r — 
l)th-order accurate approximation to q{U{xj)) by the Lipschitz continuity of qiU). Thus 
it holds that q{Uj) > e — 0{Ax'^'~~p > With the help of the concavity of q{U) and 
fl3.9p . one may know that q{Up^^^) — Eg is bounded away from zero because 


-,WENO\ 


Eg - q{U. 


h,WENO' 


, ((1 - w)Ui + w (Ui - 

^ 2At -LLF 


- en 


Fi+l - £ 


>{l-w)qiUg)+wq{Ug-j^FPl 

1 — w 

---£ — £„ > 0. 


> (1 - w)q{Uj) -Eg> 


Then, one turns to show Eg — q{Up^™^) = 0{Ax‘^'' ^). In fact, fl3.1ip implies 

QiUp-Eg = q{Up^^^)-Eg + 0{Ax) > {l-w)^-Eg + 0{Ax) > ^ (^(1 -u;)| > 0. 

Therefore 


■q{u\ 


■+,WENO^ 
j ' 


= Sr. 


g(C/+’WENO) < _ g(i7+.WENO^ ^ OiAx^'^-^] 


17 































where (13.13p has been used. 
Using the above results and 


= min 


6~t 1,6 1,6 ..il, 


gives (|3?7|) . 

The proof is completed. 


I 


3.1.2 Time discretization 


Time derivatives in the semi-discrete schemes (13.21) can be approximated by using some 
high-order strong stability preserving (SSP) method [16]. A special example considered 
here is the third order accurate SSP explicit Runge-Kutta method 


U* = U] + AtnC{U^;j), 

= 1^1 + i , 

^n+1 ^ ^ 2 ^ j)), 


with C {u-,j) = (C/) - ([/))/Ax. 


(3.14) 


In practical computations, the time stepsize selection strategy in m may be adopted 
to improve computational efficiency, and the physical-constraints-preserving flux limiter 


may also be slightly modified and implemented via enforcing directly [/** 


U' 


n+l 


G Qi and 

G Qi in the second and third stages in (13.141) . Taking the second stage as an example. 


UpWENO previous flux limiting procedure is replaced with [/±’Weno ^ _|_ 


-(u* 

4 V J 


T 


2Atr 

Ax 




3.2 Two-dimensional case 


The high-order accurate physical-constraints-preserving hnite difference WENO schemes 
presented in Subsection 13.11 can be easily extended to multidimensional RHD equations 
(II.ip . This section only presents its extension to the 2D RHD equations in the laboratory 

f 1*3^1X10 

du ^ dF,(U) , dF^W) „ 

where 


U ={D, mi, 772,2, E)'^, Fi = {Dvi, miVi + p, m 2 Vi, mi)'^, 
F2 ={Dv 2 , miV 2 , 7722X2 + P, '"22)'^. 
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Let us divide the spatial domain hi into a rectangular mesh with the cell {{x,y) \ Xj_i < 
X < < y < i/fc+i} where = (j + l)Ax and = {k + \)Ay, j,k G Z, 

and both spatial stepsizes Ax and Ay are given positive constants. 


Then a semi-discrete, (2r — l)th-order accurate, conservation finite difference scheme for 
the 2D RHD equations fl3.15l) may be derived as 


dUj,k{t) 

dt 


^l,PCP 


^l,PCP 


Ax 


+ 


-2,POP 
fc-i 


-2,PGP 


As 


=: C(U{t)-,j,k), 


(3,16) 


^ 1 PCP 2 PCP 

where the numerical flux (resp. ) is derived by using the procedure in Sub¬ 

section [3T] for each fixed k (resp. j) based on the local Lax-Firedrichs splitting and ID 
high-order WENO reconstruction with physical-constraints-preserving flux limiter. The 
time derivatives in fl3.16p may be approximated by utilizing the high-order accurate SSP 
Runge-Kutta methods, e.g. fl3.14p . 


Because of the convex decomposition 


Uj,k + AtC{U;j, k) = y (f/yfc 


2At -i,pcp 

— F -i u 
TiAx ^ 2 ' 




2At -2,PCp\ f2 




Ti f 2At -i,pcp 

+ Y 

2At -2,pcp\ 


with fi = Til (ti + T 2 ), i = 1, 2, and 


Ti = (Ax) ^max{a^.+ i J, ra = {Ay) ^max{a^.;- + i}. 


j,k 


(3.17) 


it is convenient to verify that the solutions of such resulting fully-discrete schemes belong 
to the admissible state set Q under the CFL-type condition 


At < 


w 

2(ri + T 2 )’ 


(3.18) 


where w is any positive constant less than one. In practical computations, the viscosity 
coefficients may be taken as 


Oj+l ^ max , Pi {Uj-r+l,k) , • • • , Pi {Uj+r,k)j , 

a^-fc+i = i9max |pJ;f+Y P2 {Uj,k-r+i ), • • • , P2 {Uj^k+r)^ , 

where (resp. p^’^/') is the spectral radius of the Roe matrix Ai{Uj^k, Uj+i,k) (resp. 

A 2 {Uj^k,Uj^k+i)), see [l3], approximating the Jacobian matrix Ai{U) (resp. A 2 {U)). 

The above high-order accurate finite difference schemes can also be extended to the ax- 
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isymmetric RHD equations in cylindrical coordinates (r, z) 


dU dFi{U) dF2{U) 

dt dr dz 


S{U,r), 


(3.19) 


where the flux Fi is the same as one in fl3.15p . z = 1, 2, r > 0, and the source term 


S{U,r) 


— (Dvi,miVi,m2Vi,mi)'^. 
r 


Similarly, when the computational domain hi in cylindrical coordinates (r, z) is divided 
into a uniform mesh with the rectangular cell {(r, z)\ rj_i < r < Tj+i, Zi^_i < z < 
where i = (j — |) Ar j G Z+and = (fc + |) Az, j, A; G Z, and Ar and Az are 
spatial stepsizes in r- and 2 ;-directions, respectively, the extension of the scheme fl3.16p to 
the system (I3.19p is 

-1,PCP ixl,PCP Cx2,PCP -2,POP 

dUj^k{t) _ ~ ^j+^,k ^j,k-l ~ ^3,k+^ , , . . 

~ A;^ Ai + 

=: C{U(t)-3, k) + 5([/,,fc(t), r,), (3.20) 

where the numerical fluxes ^j,k+\ ^^e same as those used in (13.1611 . 

The question is whether the scheme fl3.20p is still physical-constraints-preserving? Since 
the term Uj^k + At{C{U-,j, k) + S{Uj^k,^j)) may be decomposed into 

(1 - P)(Uj,k + k)) +P{u,,k + yS(C/,-fc,r,)), 

for any f5 G (0,1), it is sufficient to ensure that Uj^k + j,k) G Q and Uj^k + 

^S{Uj^k,rj) G G for each j,k when Uj^k ^ G for all j. A;. The first part is true if the 
condition fl3.18p is replaced with 


At < 


(1 -/3)w 

2(ri + r 2 )’ 


while the second part may be ensured if 


(3.21) 


At < BAs, As := min 


_ i^rgiUj^k) _ 

{p{U,,k) + q{U,,k))\vi{U,,k)\ 


(3.22) 


where = {{j,k) \i,k G Z,Vi{Uj^k) > 0}. The readers are referred to the following 

lemma or the similar discussion in |6l]. Combining fl3.2ip with fl3.22p . an optimal value 
of (B is chosen as u)/(m) + 2As(ri + r 2 )) such that 


(1 -/^)w 

2(ti + T 2 ) 


IBAs. 
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Lemma 3.2 IfUEQ, then U + AtS{U,r) G Q under 

_ ViAt ^ q{U) 
r ~ p + q{U)' 

Proof Assumption that U E Q implies q{U) > 0, D > 0, and p > 0. Thus if ^ < 1, then 
D{1 — ^) > 0. Hence, when .^ < 1, to ensure 

U + AtS{U,r) = {D{1 - - {E + e G, 


it is sufficient to have 

E-{E + p)^> ^{D{1 - i)f + (mi(l - 0)2 + (m 2 (l - 0)" = (1 - 0 


that is 
Therefore, if 


q{U)>i{p + q{U)). 


p + q{U) 

then U + AtS{U,r) G G- The proof is completed. 
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4 Numerical experiments 


This section conducts some numerical experiments on several ultra-relativistic RHD prob¬ 
lems with large Lorentz factor, or strong discontinuities, or low rest-mass density or pres¬ 
sure etc. to verify the accuracy, robustness and effectiveness of the proposed high-order ac¬ 
curate physical-constraints-preserving hnite difference WENO schemes. It is worth stress¬ 
ing that those ultra-relativistic RHD problems seriously challenge the numerical schemes. 
To limit the length of the paper, this section only presents the numerical results obtained 
by our hfth- and ninth-order accurate schemes with the third-order accurate Runge-Kutta 
time discretization (I3.14p . and for convenience, abbreviate them as “PCPFDWEN05” and 
“PCPFDWEN09”, respectively. Unless otherwise stated, all the computations are restricted 
to the equation of state fl2.1l) with the adiabatic index T = 5/3, and the parameter w in 
(1^ or (l3T8il is taken as 0.45 for PCPFDWEN05 and 0.4 for PCPFDWEN09. 

Example 4.1 (ID Smooth problem) This test is used to check the accuracy of our 
schemes, and similar to but more ultra than the one simulated in [SS]. The initial data 
for the ID RHD equations fl3.ll) are taken as 

U(a:, 0) = (^1-F 0.99999 sin(a;), 0.99, 0.005^ , a:G[0,27r), 

and thus the exact solution can be given as follows 

UPo:, t) = (^1-F 0.99999 sin(a: — 0.99f), 0.99, 0.005^ , xG[0,27r), t > 0. 
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It describes a RHD sine wave propagating periodically and quickly in the interval [0, 27r). 

The computational domain is divided into iVj uniform cells, i = 1,2, ••• where i is 
taken as 6 for PCPFDWEN05 and 7 for PCPFDWEN09. Here the periodic boundary conditions 
are specihed at the end points x = 0 and 27r. The time stepsize is taken as At = (0.5Aa:)3 
for PCPFDWEN05 and (0.5Ax)3 for PCPFDWEN09 in order to realize high-order accuracy in 
time in the present case. 

Tables 14.11 and 14.21 list and /°°-errors at t = 0.01 and corresponding orders obtained 
by using PCPFDWEN05 and PCPFDWEN09, respectively, where the order is calculated by 
— ln(errorj/errorj_|_i)/ln(Aj/Aj_|_i), and error* denotes the error estimated on the mesh of 
Ni uniform cell. For comparison, the errors and convergence rates are listed there for 
corresponding hnite difference WENO schemes without physical-constraints-preserving 
limiter. The results show that the theoretical order may be obtained by both PCPFDWEN05 
and PCPFDWEN09 and the physical-constraints-preserving limiter does destroy the accuracy. 

Table 4.1 


Example 14.11 Numerical errors and orders in d- and /°°-norms at t = 0.01 for PCPFDWEND5 and 
corresponding WEN05 without physical-constraints-preserving limiter. 


Ni 


WEN05 



PCPFDWEN05 


P error 

order 

l°° error 

order 

P error 

d order 

error 

order 

8 

1.8713e-3 

- 

4.4614e-4 

- 

1.8713e-3 

~ 

4.4614e-4 

- 

16 

6.7642e-5 

4.79 

1.5495e-5 

4.85 

6.7642e-5 

4.79 

1.5495e-5 

4.85 

32 

1.8277e-6 

5.21 

5.1420e-7 

4.91 

1.8277e-6 

5.21 

5.1420e-7 

4.91 

64 

5.1951e-8 

5.14 

1.6019e-8 

5.00 

5.1951e-8 

5.14 

1.6019e-8 

5.00 

128 

1.5403e-9 

5.08 

4.9554e-10 

5.01 

1.5403e-9 

5.08 

4.9554e-10 

5.01 

256 

4.6747e-ll 

5.04 

1.5215e-ll 

5.03 

4.6746e-ll 

5.04 

1.5102e-ll 

5.04 


Table 4.2 

Same as Table ITTl except for PCPFDWEN09. 


Ni 

WEN09 

PCPFDWEND9 

d error 

d order 

error 

order 

d error 

d order 

error 

order 

8 

1.2614e-4 

- 

3.0905e-5 

- 

1.2614e-4 

- 

3.0905e-5 

- 

16 

2.2845e-7 

9.11 

8.5647e-8 

8.50 

2.2845e-7 

9.11 

8.5647e-8 

8.50 

24 

5.0564e-9 

9.40 

2.3436e-9 

8.88 

5.0564e-9 

9.40 

2.3436e-9 

8.88 

32 

3.4424e-10 

9.34 

1.7915e-10 

8.94 

3.4422e-10 

9.34 

1.7915e-10 

8.94 

40 

4.3114e-ll 

9.31 

2.4253e-ll 

8.96 

4.3155e-ll 

9.31 

2.4253e-ll 

8.96 

48 

8.1007e-12 

9.17 

5.0622e-12 

8.59 

7.9810e-12 

9.26 

4.7192e-12 

8.98 

56 

1.9977e-12 

9.08 

1.1805e-12 

9.44 

1.9005e-12 

9.31 

1.1804e-12 

8.99 


Example 4.2 (ID Riemann problem) The second test is a Riemann problem (RP) 
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for the ID RHD equations fl3.1l) with initial data 


F(x,0) 


(1,0,10^)^, a:<0.5, 

(1,0,10-®)^, a:>0.5. 


(4.1) 


The initial discontinuity will evolve as a strong left-moving rarefaction wave, a quickly 
right-moving contact discontinuity, and a quickly right-moving shock wave. The flow pat¬ 
tern is similar to Example 4.2 of [58], but more extreme and difficult because of the 
appearance of the ultra-relativistic region. In the present case, the speeds of the contact 
discontinuity and the shock wave (about 0.986956 and 0.9963757 respectively) are very 
close to the speed of light. 



Fig. 4.1. Example 14.21 The density p and its close-up, the velocity vi, and the pressure p at 
t = 0.45 obtained by using PCPFDWEND5 (“*”) and PCPFDWEND9 (“o”) with 800 uniform cells. 

Fig. 14.11 displays the numerical results at t = 0.45 obtained by using PCPFDWEN05 (“*”) 
and PCPFDWEN09 (“o”) with 800 uniform cells within the domain [0,1], where the solid line 
denotes the exact solution. It can be seen that PCPFDWEN09 exhibits better resolution than 
PCPFDWEN05, and they can well capture the wave conhguration except for the extremely 
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narrow region between the contact discontinuity and the shock wave. The main reason is 
that the region between the shock wave and the contact discontinuity is extremely narrow 
(its width at t = 0.45 is about 0.00424) so that it can not be well resolved with 800 uniform 
cells. At the resolution of 800 uniform cells, the maximal densities for PCPFDWEN05 and 
PCPFDWEN09 within the above narrow region are about 58.7% and 74.4% of the analytic 
value, respectively. 


Example 4.3 (Blast wave interaction) This is an initial-boundary-value problem for 
the ID RHD equations fl3.ip and has been studied in [331158] . The same initial setup is 
considered here. The adiabatic index T is taken as 1.4, the initial data are taken as follows 


E(x,0) 


'( 1 , 0 , 1000 )^, 
< ( 1 , 0 , 0 . 01 )^, 
.( 1 , 0 , 100 )^, 


0 < X < 0.1, 

0.1 < X < 0.9, 
0.9 < X < 1, 


(4.2) 


and outflow boundary conditions are specihed at the two ends of the unit interval [0,1]. 
This is also a very severe test since it contains the most challenging one-dimensional rel¬ 
ativistic wave conhguration, e.g., strong relativistic shock waves, and interaction between 
blast waves in a narrow region, etc. 


Fig. 14.21 gives close-up of the solutions at f = 0.43 obtained by using PCPFDWEN05 (“*”) 
and PCPFDWEN09 (“o”) with 4000 uniform cells within the domain [0,1]. It is found that 
the solutions at f = 0.43 within the interval [0.5, 0.53] consists two shock waves and 
two contact discontinuities since both initial discontinuities evolve and both blast waves 
collide each other; and compared to the GRP scheme in [58| , both schemes can well resolve 
those discontinuities and clearly capture the complex relativistic wave conhguration, but 
PCPFDWEN09 exhibits better resolution than PCPFDWEN05 except for slight overshoot and 
undershoot of the rest-mass density between the left shock and the contact discontinuity. 
The overshoot and undershoot may be suppressed by using the monotonicity-preserving 
limiter in [2], see Fig. 14.31 

Example 4.4 (Shock heating problem) The last ID test is to solve the shock heating 
problem, see [3], by using PCPFDWEN05 and PCPFDWEN09. The computational domain [0,1] 
with a rehecting boundary at the right end is initially hlled with a cold gas (the specihc 
internal energy is taken as 0.0001). The gas has an unit rest-mass density, the adiabatic 
index F of 4/3, and the velocity xq of 1 ~ 10“^°, When the initial gas moves toward to the 
reflecting boundary, the gas is compressed and heated as the kinetic energy is converted 
into the internal energy. After then, a reflected strong shock wave is formed and propagate 
to the left with the speed Xg = (F — l)lFo|xo|/(lFo + 1), where Wq = (1 — Xg)”^/^ is about 
70710.675. Behind the reflected shock wave, the gas is at rest and has a specihc internal 
energy of Wq — 1 due to the energy conservation across the shock wave. The compression 
ratio a across the relativistic shock wave 

^ - 1) ~ 282845.7, 

grows linearly with the Lorentz factor Wq and towards to inhnite as the inhowing gas 
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120 


(a) p 


(b) vi 


0.505 0.51 0.515 0.52 0.525 0.53 



(c) p 


(d) e 


Fig. 4.2. Example 14.31 Close-up of the numerical solutions at t = 0.43 obtained by using 
PCPFDWEND5 (“*”) and PCPFDWEND9 (“o”) with 4000 uniform cells. The solid lines denote the 
exact solutions. 


velocity approaches to speed of light. It is worth noting that the compression ratio across 
the shock wave in the non-relativistic case is always bounded by (F + l)/(r — 1 ). 


Fig. 14.41 displays the numerical solutions at t = 2 obtained by using PCPFDWEN05 (“*”) 
and PCPFDWEN09 (“o”) with 200 uniform cells. It shows that both schemes exhibit good 
robustness for this ultra-relativistic problem and good resolution for the strong shock 
wave, even though there exist slight oscillations in the rest-mass density and the internal 
energy behind the shock wave as well as well-known wall-heating phenomenon near the 
reflecting boundary x = 1. Similar to Example 14.31 those small oscillations can also be 
efficiently alleviated by using the monotonicity-preserving limiter [2]. To save the paper 
space, the results are not presented here. 


Besides the flux limiter in Section 13.1.11 we have also tried to extend the parametrized 
flux limiter [5^291154] to the ID RHD equations fl3.ip . the numerical results show that 
the parametrized flux limiter is less restrictive on the CFL number in preserving high 
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(c) p (d) e 

Fig. 4.3. Same as Fig. 14.21 except for PCPFDWEND9 (“o”) with a monotonicity-preserving limiter. 

order accuracy and has slightly better resolution for Example 14.21 than the flux limiter in 
Section 13.1.11 


Example 4.5 (2D Riemann problem) The non-relativistic 2D RPs are theoretically 
studied for the first time in [HB]. After then, the 2D RPs become benchmark tests for 
verifying the accuracy and resolution of numerical schemes, see 


Initial data of two RPs of 2D RHD equations fl3.15p considered here comprise four different 
constant states in the unit square D = [0,1] x [0,1], while initial discontinuities parallel 
to both coordinate axes respectively. 


The initial data of the hrst RP are 


V{x,y,0) 


'( 0 . 1 , 0 , 0 , 0 . 01 )'^, 
(0.1,0.99,0,1)^, 
(0.5, 0,0,1)^, 
^(0.1,0,0.99,1)^ 


X > 0.5, y > 0.5, 
X < 0.5, y > 0.5, 
X < 0.5, y < 0.5, 
X > 0.5, y < 0.5, 
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Fig. 4.4. Example 14.41 The density p, the velocity vi, the pressure p, and the internal energy 
e at t = 2 obtained by using PCPFDWEN05 (“*”) and PCPFDWEN09 (“o”) with 200 uniform cells. 
The solid lines denote the exact solutions. 

where both the left and bottom discontinuities are contact discontinuities with a jump in 
the transverse velocity, while both the right and top discontinuities are not simple waves. 
Note that this test is different from the case in |59j . 


Fig. 14.51 gives the contours of the density logarithm In p at time t = 0.4 obtained by 
using PCPFDWEN05 and PCPFDWEN09 with several different mesh resolutions. It is found 
that four initial discontinuities interact each other and form two reflected curved shock 
waves, an elongated jet-like spike, which is approximately between two points (0.7,0.7) 
and (0.9,0.9) on the diagonal x = y when t = 0.4, and a complex mushroom structure 
starting from the point (0.5,0.5) respectively and expanding to the bottom-left region; 
both PCPFDWEN05 and PCPFDWEN09 well capture these complex wave configuration, and it 
is obvious that with the same mesh of 400 x 400 uniform cells, PCPFDWEN09 gets better 
resolution of discontinuities than PCPFDWEN05, and the solutions obtained by PCPFDWEN09 
with the mesh of 800 x 800 uniform cells are comparable to those obtained by PCPFDWEN05 
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with a finer mesh of 1200 x 1200 uniform cells. For a further comparison, the rest-mass 
densities are plotted along the line y = x and y = I, see Fig. 14.61 Those plots validate the 
above observation. 



(a) PCPFDWEND5 with 400 x 400 uniform cells (b) PCPFDWEN09 with 400 x 400 uniform cells 



(c) PCPFDWEND5 with 1200 x 1200 uniform cells (d) PCPFDWEN09 with 800 x 800 uniform cells 


Fig. 4.5. The first 2D RP in Example 14.51 The contours of the density logarithm Inp at t = 0.4 
within the domain [0,1] x [0,1] obtained by using PCPFDWEN05 and PCPFDWEN09 (25 equally 
spaced contour lines from —6 to 1.9). 
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Fig. 4.6. Same as Fig. 14.51 except for p along the line y = x within the scaled interval [0,1] (left) 
and p along the line y = 1 within the closed interval [0.75,0.95] (right). 

The initial data of the second 2D RP are 


F(x,y,0) 


'( 0 . 1 , 0 , 0 , 20 )^ 

(0.00414329639576, 0.9946418833556542, 0,0.05)^, 
(0.01,0,0,0.05)^, 

^ (0.00414329639576, 0, 0.9946418833556542, 0.05)^, 


X > 0.5, y > 0.5, 
X < 0.5, y > 0.5, 
X < 0.5, y < 0.5, 
X > 0.5, y < 0.5, 


where both the left and bottom discontinuities are contact discontinuities while both the 
top and right are shock waves with the speed of —0.66525606186639. As the time increases, 
the maximal value of the fluid velocity may be very close to the speed of light. 

Fig. 14.71 displays the contours of the density logarithm Inp at time t = 0.4 within the 
unit square [0,1] x [0,1] obtained by using PCPFDWEN05 and PCPFDWEN09 with the mesh 
of 400 X 400 uniform cells. It is seen that the interaction of four initial discontinuities 
leads to the distortion of both initial shock wave and the formation of a “mushroom 
cloud” starting from the point (0.5, 0.5) and expanding to the left bottom region, and the 
present schemes have good performance and strong robustness in such ultra-relativistic 
flow simulation, in which the fluid velocity reachs 0.9998458 locally, or the Lorentz factor 
may be not lesser than 56.95. Plots of Inp along the line y = x in Fig. 14.81 further show 
that PCPFDWEN09 captures the structure better than PCPFDWEN05. 

Example 4.6 (Forward facing step problem) The forward facing step problem was 
hrst introduced by Emery [T2], has been widely used to test the non-relativistic hydrody¬ 
namic codes, e.g. |1], and extended to the ideal relativistic fluid, see [5DEI] - 

The same setup as in ISD is used here. The wind tunnel is located in the domain [0, 3] x [0,1] 
and contains a forward facing step with a height of 0.2, which starts from x = 0.6 and 
continues along the length of the tunnel. Initially it is hlled with an ideal gas with the 
density of 1.4, the velocity of 0.999, the Mach number of 3, and the adiabatic index of 
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Fig. 4.7. The second 2D RP in Example ld.St The contours of the density logarithm In p at t = 0.4 
within the domain [0,1] x [0,1] obtained by using PCPFDWEN05 and PCPFDWEN09 with 400 x 400 
uniform cells (25 equally spaced contour lines from —8 to —2.3). 



Fig. 4.8. Same as Fig. 14.71 except for Inp along the line y = x within the scaled interval [0,1]. 

1.4. The reflective boundary conditions are specified on the walls of the tunnel, while the 
inflow and outflow boundary conditions are specified at two ends of the tunnel (x = 0 
and 3). Due to the forward facing step, the bow shock wave is formed and then a Mach 
reflection happens at the top wall of the tunnel. Later, the reflected shock wave incidents 
to the step and then the regular reflection of the shock wave is caused and generates a 
second reflected curved shock wave. Moreover, the step corner (0.6,0.2) is the center of 
a rarefaction fan and hence a singular point of the flow. Although the time-evolution of 
the flow is similar to the non-relativistic case, see e.g. [1], the present test is much more 
difficult than the non-relativistic case because the bow shock wave has a much fast speed 
and there exist the large jumps in the density and the pressure and the ultra-relativistic 
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regime where the Lorentz factor hh S> 1. 



Fig. 4.9. Example 14.61 The contours of the rest-mass density logarithm Inp at t = 4 obtained by 
using PCPFDWEN05 (top) and PCPFDWEND9 (bottom) with 300 x 100 uniform cells for the domain 
[0, 3] X [0,1] (25 equally spaced contour lines from -0.86 to 4.64). 

Fig. 14.91 gives the contours of the rest-mass density logarithm Inp at f = 4 obtained by 
using the proposed PCPFDWEN05 and PCPFDWEN09 with 300 x 100 uniform cells for the 
domain [0, 3] x [0,1]. Numerical results also show that across the Mach stem, the relative 
jumps Ap = \pr/Pl — 1| and Ap = \pr/pl — 1| are about 61.33 and 5223.99, respectively, 
where the subscripts L and R denote corresponding left and right states. Comparing our 
results with those in [HI], it can be seen that the flow structures are well obtained and 
the discontinuities, e.g. the bow and reflected shock waves and the Mach stem, are well 
captured with high resolution by using PCPFDWEN05 and PCPFDWEN09 without any artihcal 
entropy £x around the step corner. 

Example 4.7 (Axisymmetric relativistic jets) The last 2D example is to simulate 
two high-speed relativistic jet flows by solving the axisymmetric RHD equations fl3.19p . 
The jet flows are ubiquitous in extragalactic radio sources associated with active galactic 
nuclei and the most compelling case for a special relativistic phenomenon. Since there are 
the ultra-relativistic region, strong relativistic shock wave and shear flow, and interface 
instabilities etc. in the high-speed jet flow, simulating successfully such jet flow can be a 
real challenge, see e.g. [361I11II35II26] . 

The hrst test is a pressure-matched hot Al model, see [35]. In this model, the classical 
beam Mach number Mi, is near the minimum Mach number M™" = vi,/\/T — 1, the beam 
is moving at speed Vb = 0.99, and relativistic effects from large beam internal energies 
are important and comparable to the effects from fluid velocity near the speed of light. 
Initially, the computational domain [0,7] x [0, 50] in the (r, z) plane is hlled with a static 
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Fig. 4.10. The hot Al model in Example 14.7t Schlieren images of the rest-mass density logarithm 
ln/9 within the symmetrical domain [—7,7] x [0,50] at several different times obtained by using 
PCPFDWEN05 with 280 x 2000 uniform cells. 
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Fig. 4.11. Same as Fig. 14.101 except for PCPFDWEND9. 

uniform medium with unit rest-mass density with the adiabatic index of 4/3. A light 
relativistic jet is injected in the ^-direction through the inlet part (r < 1) of the bottom 
boundary {z = 0) with a density of 0.01, a pressure equal to the ambient pressure, and a 
speed of Vb or a Lorentz factor of 7.09. The relativistic Mach number Mr := MbW/Ws is 
about 9.97, where IT/ = 1 /y/l — is the Lorentz factor associated with the local sound 
speed and Mb = Vb/cs is equal to 1.72. The symmetrical boundary condition is specihed 
at r = 0, the hxed inflow beam condition is specihed on the nozzle {z = 0,r < 1}, while 
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outflow boundary conditions are on other boundaries. 


Figs. 14.101 and 14.111 display the schlieren images of the rest-mass density logarithm In p 
within the domain [—7,7] x [0,50] at t = 10,20,30,40,50 and 60 obtained by using 
PCPFDWEN05 and PCPFDWEN09 with 280 x 2000 uniform cells respectively. Compared them 
to those in [55], it is found that the time evolution of a light, relativistic jet with large 
internal energy is well simulated by our schemes, and the Mach shock wave at the jet head 
is correctly and well captured during the whole simulation. Although the internal structure 
is almost completely lacked because of the pressure equilibrium between the beam and 
its surroundings, and the beam/cocoon interface of the hot jets is very stable against the 
growth of pinch instabilities that would evolve into internal shock waves, the proposed 
schemes still clearly resolve the beam/cocoon interface and the Kelvin-Helmholtz type 
instability at the beam/cocoon interface. 



Fig. 4.12. The highly supersonic jet test of Example l4.7l Schlieren images of the rest-mass density 
logarithm Inp at t = 100 obtained by PCPFDWEN05 (left) and PCPFDWEM09 (right) on the mesh of 
384 X 1152 uniform cells. 

The second test is the pressure-matched highly supersonic C2 jet model {Mb S> 
r = 5/3, and vt = 0.99), see [351I6T] . Highly supersonic jets are also refered to cold models 
and the relativistic effects from large beam speeds dominate in C2 jet model so that there 
exist important differences between hot and cold relativistic jets. 

The same setup as in [6T] is used here within the computational domain [0,15] x [0,45] 
in the (r, z) plane. The initial relativistic jet has the same rest-mass density and velocity 
as those in the above hot Al model, but a higher Mach number of 6 (or corresponding 
relativistic Mach number of about 41.95) and a larger adiabatic index of 5/3. 

Fig. 14.121 shows the schlieren images of the rest-mass density logarithm Inp within the 
symmetrical domain [—15,15] x [0,45] at t = 100 obtained by using PCPFDWEN05 and 
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PCPFDWEN09 with 384 x 1152 uniform cells. In our simulations, stable cocoons can be found 
in the early stages of the evolution, but these cocoons eventually evolve into large vortices 
producing turbulent structures, and the morphology and dynamics of the relativistic jets 
in Fig. 14.121 agree well with those obtained by using an adaptive mesh refinement RHD 
code in [6T]. It is worth noting that the above-observed gross morphological features 
already found in non-relativistic calculations. 


5 Conclusions 


The paper developed (2r — l)th-order accurate finite difference WENO schemes for special 
RHD equations and proved that their solutions satisfied physical properties: the positivity 
of the rest-mass density and the pressure and the bounds of the velocity. The key con¬ 
tributions were that the convexity and some mathematical properties of the admissible 
state set Q were proved and the concave function q{U) with respect to the conservative 
vector U was discovered. The schemes were built on the local Lax-Friedrich splitting, the 
WENO reconstruction, the physical-constraints-preserving flux limiter, and the high-order 
strong stability preserving time discretization. They were considered as formal extensions 
of positivity-preserving finite difference WENO schemes for the non-relativistic Euler 
equations [21]. However, developing physical-constraints-preserving methods for the RHD 
system became much more difficult than the non-relativistic case because of the strongly 
coupling between the RHD equations, no explicit expressions of the primitive variables 
V and the flux vector Fi in terms of the conservative vector U, and one more physi¬ 
cal constraint for the fluid velocity in addition to the positivity of the rest-mass density 
and the pressure, the non concavity of p{U), which was a key ingredient to enforce the 
positivity-preserving property for the non-relativistic Euler equations. 

Several numerical examples demonstrated accuracy, robustness, and effectiveness of the 
proposed physical-constraints-preserving schemes in solving relativistic problems with 
large Lorentz factor, or strong discontinuities, or low rest-mass density or pressure etc., 
which involved a ID smooth problem, a ID Riemann problem, a ID blast wave interac¬ 
tion problem, and a ID shock heating problem, two 2D Riemann problems, and a forward 
facing step problem and two axisymmetric jet flows in two dimensions. 


Since the present physical-constraints-preserving limiting procedure was independent on 
the reconstruction in space, other high-order accurate reconstruction or interpolation 
could also be used to replace the WENO reconstruction. Moreover, it was possible that the 
analysis results and the limiting procedure etc. could be extended to develop high-order 
accurate physical-constraints-preserving finite volume schemes for the RHD equations. 
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